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While quantum computers are capable of simulating many quantum systems efficiently, the simulation al- 
gorithms must begin with the preparation of an appropriate initial state. We present a method for generating 
physically relevant quantum states on a lattice in real space. In particular, the present algorithm is able to prepare 
general pure and mixed many-particle states of any number of particles. It relies on a procedure for converting 
from a second-quantized state to its first-quantized counterpart. The algorithm is efficient in that it operates in 
time that is polynomial in all the essential descriptors of the system, such the number of particles, the resolution 
of the lattice, and the inverse of the maximum final error. This scaling holds under the assumption that the 
wavefunction to be prepared is bounded or its indefinite integral known and that the Fock operator of the system 
is efficiently simulatable. 



I. INTRODUCTION 

Simulating quantum systems on a conventional computer 
requires resources that generally scale exponentially with the 
size of the system. Feynman proposed to solve this problem 
using a quantum machine that would be able to mimic the 
properties of the quantum system [1]. Subsequently, it has 
been demonstrated that quantum computers would be able to 
simulate the time-dependent Schrodinger equation for many 
systems of interest using resources that scale polynomially 
with the size of the system H d S II H EE S H • 

However, all 

such simulations require the preparation of an appropriate ini- 
tial state, which must be preparable to within a chosen error. 

In this work, we focus on the preparation of states on a gate- 
model quantum computer. Our techniques can therefore com- 
plement those developed for the preparation of states in other 
models of quantum computation, such as adiabatic quantum 
computing @,[T3]. 

In general, we will call a state on n qubits "efficiently 
preparable" if it can be prepared, to within error e, using 
poly(n, e _1 ) elementary (one- and two-qubit) quantum gates. 
Unfortunately, the efficiently preparable states form only a 
small subset of all quantum states. This is because a gen- 
eral state on n qubits contains 2" amplitudes, and therefore 
one needs 0(2") gates to prepare it [11]. Indeed, state- 
preparation algorithms are known that almost reach this lower 
bound lUlHai. 

In this work, we show that if wavefunctions are represented 
on a grid in real space, then most quantum states of physi- 
cal interest are efficiently preparable. This is of interest be- 
cause efficient, grid-based simulation algorithms are known 
for physically realistic systems JH 0, HI 

Our work extends that of Zalka, who, in introducing real- 
space quantum simulation ||4j], also provided the first state 
preparation algorithm. However, his procedure is able to 
prepare only states of single particles or uncorrected many- 
particle systems. We show how to use Zalka's single-particle 
wavefunctions as building blocks, permitting the preparation 
of general superpositions and mixed states of an arbitrary 
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number of particles. Our approach is motivated by electronic- 
structure theory, in that we we choose particularly convenient 
single-particle bases in which to expand more complicated 
states. We use the single-particle eigenstates to form Slater de- 
terminants (configurations), superpositions of which are used 
to express general many-particle states. 

Our scheme is essentially a method for translating states in 
second quantization to the corresponding states in first quan- 
tization. This has two advantages. First, many useful states 
that might be needed in first quantization are easily prepared 
in second quantization flU . In particular, we can prepare 
eigenstates of operators if our scheme is combined with full 
configuration interaction (FCI) [15], an exact diagonalization 
method. FCI is classically an exponentially hard problem 
due to the exponential growth of the number of configura- 
tions with system size, but it can be computed on a quantum 
computer in polynomial time [6]. The quantum FCI operates 
in second quantization, and can compute, for example, the 
ground state wavefunction of a molecular system. The sec- 
ond benefit of our method is that it is often easier to simulate 
time-evolution in real space than in Fock space. For instance, 
every simulation in second quantization would require a sep- 
arate set of basis-set-dependent operators and there might be 
some processes, such as ionization, which could not be ade- 
quately described using a small, localized basis set. In first- 
quantization, however, all problems of chemical interest can 
be efficiently simulated by direct simulation of the molecular 
Hamiltonian in real space |8||. 

This paper is organized as follows. We first consider the 
preparation of many-particle states in which all the particles 
are the same. There are three steps: the preparation of single- 
particle eigenstates in a chosen basis, the preparation of many- 
particle configurations, and finally the preparation of superpo- 
sitions of configurations. We discuss the preparation of mixed 
states, after which we turn to systems with many different 
types of particles. Again, we consider the preparation of con- 
figurations, their superpositions, and mixed states. We close 
by showing that the algorithm is efficient in that its run-time 
is polynomial in the size of the system, the number of qubits 
used to encode the wavefunction, and the inverse of the maxi- 
mum allowed error. 
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II. ONE TYPE OF PARTICLE 

Our algorithm translates from second to first quantization, 
and therefore depends on the basis which is chosen for the 
representation of the second-quantized states. We require a 
finite orthonormal basis of functions {<pi, . . . , 4>ai}, which are 
the eigenstates of a known operator F on an A/-dimensional, 
single-particle Hilbert space. In our analogy with electronic- 
structure theory, F would be the Fock operator for a single 
particle lll5ll . and indeed we expect that the algorithm would 
be at its most useful if F is chosen as the Fock operator of 
an actual system. Although the form of F can be arbitrary, 
subject to a few restrictions below, we will take advantage of 
the analogy and refer to F as the Fock operator. We will also 
say, for example, that two eigenstates of F are degenerate if 
their energies (the corresponding eigenvalues) are the same. 

To ensure that the overall state-preparation algorithm scales 
efficiently, we require that F can be efficiently simulated 
on a quantum computer, i.e., that the simulation time scales 
polynomially with the size of the system. More precisely, if 
there are m particles occupying the M orbitals and the sim- 
ulation is done on a grid of 2 l sites (see below), then, for 
any t and any e > 0, there exist a unitary U, composed of 
poly(m, M, I, t, e^ 1 ) elementary (one- and two-qubit) quan- 



tum gates, such that 



< e. Intuitively, this means 



U - e~ tFt 

that given an initial state, the final state generated by the ac- 
tion of F for time t can be calculated with reasonable effort 
and reasonable error. 

Several classes of Hamiltonians are known to be efficiently 
simulatable, and together they ensure that most physically rel- 
evant Fock operators will also be efficiently simulatable. Very 
generally, an operator can be efficiently simulated if its matrix 
in a given basis is sparse [fToL HH [13] • in particular, this in- 
cludes Hamiltonians that are sums of local operators, each of 
which acts on only a few degrees of freedom JH Ell • In addi- 
tion, many physically realistic real-space Hamiltonians (such 
as those for chemical systems) can be efficiently simulated 



We finally note that the requirement that the basis be or- 
thonormal may exclude certain commonly used basis sets. 
Many of the usually encountered bases are appropriate, such 
as plane waves or molecular orbitals, which diagonalize 
the molecular Hartree-Fock Hamiltonian. However, non- 
orthogonal bases, such as Gaussian wavepackets or atomic or- 
bitals on more than one atomic center, are not suitable for state 
preparation using our procedure. 



A. Single-particle eigenstates 

A single-particle basis function <\> can be prepared on a grid 
by the state preparation method first proposed by Zalka [4] 
and rediscovered independently by both Grover and Rudolph 
lfl8tl and Kaye and Mosca fl9jl. The algorithm first prepares 
the absolute value of the function, followed by the addition of 
the phases. Specifically, given a register of I qubits, represent- 



ing a grid of 2 l points, and a basis state <f>(x) normalized over 
a length L, the algorithm first performs the transformation 



io> - 10> = E 



\x), 



where each integer- valued state \x) is a position on the suit- 
ably scaled grid. This state is generated from the state 
1 000...) by redistributing its amplitude I times across the 
eigenstates \x). To perform the redistribution correctly, we 
calculate the integrals 



Ii,k — 



f^\^{x)\Hx 

r, ( i +2) ^)iW 



(i) 



for k — 0, . . . , 2* — 2 and i — 1, . . . , I, The fraction J,-^ is 
simply the probability that a particle in the (fc + l)th subdi- 
vision of size L/2 1 is also in its left half. If the denominator 
in Ii t k is zero, there is no amplitude to redistribute, so we can 
skip this step. The first split is realized by performing a rota- 
tion on the first qubit by arccos(y / /i.o), corresponding to the 
transformation 



10,. 



VT^o|l, 



This splits the norm of the initial state so that the appropriate 
proportion is present on each half of the grid. The subsequent 
finer splits are carried out in superposition using controlled 
rotations on each qubit. For example, after the second iter- 
ation, the correct proportion of the norm is present in each 
quarter of the grid. After I iterations, one obtains the desired 
state. Note that adding a single qubit and the corresponding 
rotation doubles the precision of the grid. Consequently, the 
absolute value of the wavefunction can be efficiently approx- 
imated to any desired accuracy. Phases can be added where 
necessary through phase-kickback [20]. Given a procedure 



that can transform \x) 
preparation of | <j>) as 



x), we can complete the 



E 
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The same algorithm can be straightforwardly generalized 
to a three-dimensional grid, where the position eigenstates 
are in Cartesian coordinates and the corresponding three- 
dimensional integrals are used. In addition, particle spin can 
be represented using additional qubits. A particle with spin S 
requires [log 2 (2S' + 1)] qubits to store its z-projection ms- 
In particular, there is a natural mapping between the spin of 
spin-i particles and the states of a single qubit. If the Fock 
operator is spin-free, the eigenstates will have separable spa- 
tial and spin degrees of freedom, making the complete single- 
particle state \4>)\ms). Preparing the spin part of this wave- 
function is relatively easy, for it suffices to initialize the spin 
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register in an integer state. If, however, the eigenstate has 
correlation between the spatial and spin degrees of freedom, 
it can be prepared using the techniques in Sees. Ill Dl and Hill 
That is, we treat the particle as if it were a composite system — 
composed of a spinless, spatial part and a spin — and prepare 
its eigenstate using the techniques below. In what follows, we 
will assume that our particles are fermions and we will note 
where the algorithm would need to be modified for bosons. 



B. Computational complexity of integration 

The preceding method for preparing single-particle states 
requires the evaluation of integrals (Q~|). Since this must be per- 
formed in superposition, the integrals must be computed on 
the quantum computer: precomputing them classically would 
require an exponentially large look-up table. Consequently, 
the computational complexity of the state preparation proce- 
dure will depend on the the cost of computing the integrals 

EM. 

An integration procedure will, given a function <fi : V C 
R d — > E (where V is a bounded region), supply an estimate I 
of the integral I = J 0(x)e?x such that \I — I\ < £/ with a 
certain fixed probability 5 (we'll call this condition the (sj, S) 
absolute error). 

Integrals can be evaluated either analytically or numeri- 
cally. If the indefinite integrals of the basis functions are 
known, the definite integrals over any box on the Cartesian 
grid can be computed. The values of the indefinite integrals 
themselves can usually be computed efficiently (i.e., with 
polynomial cost in the desired accuracy) because they usu- 
ally contain simple mathematical functions. In particular, the 
time it takes to retrieve n digits of any elementary function 
is a polynomial in n i2lll . and likewise for compositions of 
elementary functions. 

If the indefinite integrals are either unknown or impractical 
to compute, numerical techniques can be used. In particular, 
any classical numerical technique can, in principle, be imple- 
mented on a quantum computer. For example, computing I 
by Monte Carlo requires, in the worst case j22ll . 



(S-^l-f/^V/e? 



samples of (f> for an (£/,£) absolute error, 
an estimate of the variance of over V and $(z) = 
(27r) -1 / 2 J* exp(— u 2 /2)du is the standard normal cumu- 
lative distribution function. In particular, if <p is bounded so 
that </>l < 4>(x) < (jm for all x £ V, the number of required 
samples is limited 112211 to 



Furthermore, it is known that quantum computers are able 
to offer a quadratic speed-up over conventional probabilistic 
methods of integral evaluation. Quantum integration tech- 
niques lf24ll25ll rely on amplitude amplification Iu6ll to achieve 
a computational complexity of 0{e^ 1 ). This has been proven 
optimal by Nayak and Wu Il27[|28ll . These techniques have the 
same general applicability as classical Monte Carlo, and will 
likewise succeed for any bounded function. Furthermore, the 
state preparation scheme of Soklakov and Schack [29], which 
relies on amplitude amplification, also succeeds in 0(ej 1 ) 
time. 

The preceding assumes that the function that we seek to 
prepare does not depend substantially on the grid spacing. We 
would expect that of realistic wavefunctions, assuming that 
the grid spacing is smaller than the smallest wavelength of the 
system. A useful exception are Kronecker (5-functions, de- 
fined on a grid of 2' points as </>(x) = 2'/ 2 <5 X Xo , where Xo is a 
constant vector. The variance of (5-functions grows exponen- 
tially in I, and therefore they cannot be integrated efficiently 
by Monte Carlo or prepared efficiently using the method of 
Soklakov and Schack II29I1 . However, they can still be pre- 
pared efficiently using our techniques because their indefinite 
integral, the Heaviside function, can be easily computed in 
time independent of 

It remains to be shown that an error in the evaluated integral 
translates to a comparable error in the prepared function. If the 

integrals ([TJ have a maximum error si, that is, 1^ — In- < 



Ei, and the error in the final prepared state </>y is = 1 — 
l^(f> <fij , we find that < lej/2. In the case I — 1, \<p) — 
-\/J|0> + x/T - 7|1) and 4>) = y/l \0) + Then, 



= V7|C 

assuming < I < 1, which is necessary for 
acceptable state, 



to be an 



£0 = 1-VJJ- V(i-/)(i--0 



< 1 - 
£1 

~ 2 



where we have assumed that ej -C 1. For larger I, a sim- 



ilar analysis applies qubit-wise: one finds that ((f>cf)) > 

(1 — £/) ' > 1 — lej/2 (the last inequality holds for all sj 
if I > 2), whence < lei/2. That is to say, the error in the 
prepared state grows only polynomially with the error in the 
evaluated integral, a fact that we will use later on to establish 
the computational cost of the state preparation algorithm. 



{^-\\-S/2){ih-H)l^i)' 



That is, Monte Carlo integration of any finite-variance func- 
tion requires time that scales as 0(ej 2 ). Acceptable wave- 
functions need not be continuous or even finite 112 311 (and 
hence may have infinite variance), but such examples are 
rather contrived and rarely encountered in practice (but see 
below for (5-functions). 



C. Many-particle eigenstates 

The next step is to use Zalka's algorithm to prepare multi- 
particle configurations. That is, we wish to prepare the 
position-space representation of a second-quantization state 
\n\n2 ■ ■ ■ n>M)2nd ( a F° c k eigenstate), where n, is the occu- 
pation number of the basis orbital (pi (1 < i < M). The 
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position-space representation of \nin 2 . . . nM) 2n d w ^ be a 
Slater determinant of the occupied orbitals, and it will be an 
eigenstate of the many -body Hartree-Fock Hamiltonian 



H 



(2) 



where m is the total number of particles and Fi is the Fock 
operator F acting on the particle i lfl5ll . 

We assume that the state \n1n2 ■ - ■ nM} 2n d ^ as au * ea dy 
been prepared by some previous algorithm. The M basis or- 
bitals are occupied by m particles and we let j\ , . . . , j m be the 
indices of the occupied orbitals. We therefore wish to perform 
the transformation 

I nm2 ... nM) 2 nd <8> |0 ... )i» t — > 
|nin 2 ...nAf)2ni®-7= s sK (7 )l cr ( < fe < fe ■ ■ ■ <t>j m ))ut 

= \nin 2 ■ ■ ■ n M ) 2nd ® \$) lst , (3) 

which takes the input state and prepares the appropriate first- 
quantized Slater determinant |$) lst , a superposition of all the 
permutations on the m occupied orbitals (S m is the symmet- 
ric group on m elements and sgn denotes the signature). Here 
|0 . . . )i s t contains m registers for the m first-quantized oc- 
cupied orbitals \<j>j 1 ) lst , ■ • ■ , \ <Pj m )± +■ Note that <[3j is not in 
general a reversible operation, as multiple input states would 
be mapped to the same antisymmetrized result. To ensure 
the algorithm is reversible, we additionally require fl3Qfl that 
ji < J2 < • • • < jm- The procedure can be slightly modified 
if bosons are in question: then sgn(cr) is to be omitted, and 
the ji must satisfy ji < J2 < ■ ■ ■ < jm- 

The transformation (0) is accomplished in two steps. First, 
the occupied single-particle basis orbitals are each prepared in 
a separate register, forming a Hartree product: 



\nin 2 ■ ■ ■ n M )2nd ® |0 . . . )i sf -> 
\nin 2 - - -n M )2nd ' 



list- 



The procedure can be modified in the case of bosons by count- 
ing the occupation of each orbital and preparing that many 
copies in separate registers. 

In the next step, the Hartree product is antisymmetrized, 
which produces the desired Slater determinant. To complete 
this step, we introduce an improved form of the antisym- 
metrization algorithm developed by Abrams and Lloyd IB0I1 . 
The algorithm begins with the m wavefunctions | <f>j i ) to be 
antisymmetrized in register A, and m [log 2 m~\ qubits in reg- 
ister B (where each grouping of [log 2 m~\ qubits constitutes 
a "quword") initialized to |0). Using a series of controlled 
rotations, B is converted to the state 



^ m m — 1 



[2] 



|1) 



B\m] 



i=l 



which is a superposition of m! unique states consisting of m 
quwords each, and B[i] denotes the ith quword in register B. 



Next we will transform this state into the superposition 
7= \<r{l,...,m)) B 



<J£S„ 



as follows. First let B'[l] = B[l]. Then assign to 
B'[i] the _B[z]th natural number not present in the set 
{B'[1],B'[2}, B'[i - 1]}. This leaves the quantum com- 
puter in the state 



1 



\a{l,...,m)) B . (4) 



Register B now contains a symmetrized state and this sym- 
metry can be transferred to register A by sorting B while 
performing the same swaps on the wavefunctions in A. This 
yields the symmetrized state 



1 



cr£S m 



|c(<An<A?2 ■■■<f>j m ))A® |l,2,...,m> B , 



(5) 



which is what we would keep if we were interested in prepar- 
ing bosonic states. To instead obtain an antisymmetrized state, 
we need only count the number of exchanges made in the sort, 
and reverse the sign of the wavefunction if it is odd. If we now 
eliminate the register B, A contains the desired multi-particle 
state |$)i si . 

The original algorithm, introduced by Abrams and Lloyd, 
included an additional auxilliary register C, which would then 
be used as an intermediate for the sorting of A and B. We 
eliminate this step by sorting A and B together directly. 



D. Superpositions 

We now generalize the algorithm to the preparation of su- 
perpositions of many-particle states. Given a superposition 
of second-quantization states |n i ) 2 „d = \nun 2i . . . n M i)<2nd, 
with amplitudes a>i, we wish to perform the transformation 



ai\ni) 2n d ®|0. 



list 



)2nd® 



The superposition on the left might come from a variety of 
sources. For example, an easily-prepared equal superposition 
of Fock states would result in an equal superposition of real- 
space wavefunctions. Wang et al. provide an algorithm for 
preparing general superpositions of Fock states on a quantum 
computer 11411 . Alternatively, a quantum electronic-structure 
algorithm could be used to efficiently produce a physically 
relevant superposition. For example, an FCI algorithm could 
specify the ground state of a chemical or other many-body 
system in terms of a superposition of Fock states Jg]. 

As before, we begin by applying Zalka's state preparation 
algorithm to the input state. Because this linear operation is 
carried out in superposition, it accomplishes the transforma- 
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tion 

^^a i |n i ) 2 „d^ ® |0 . . .) lst -> 

2J Q!i|ni)2«d <8> \<f>jii<f>j 2 i ■ ■ ■ <t>j m i)lst- 

i 

Note that the single-particle wavefunctions are now entangled 
with the input state. For a multi-particle eigenstate, the sit- 
uation was different because the resulting state was separa- 
ble. Hence, to separate the first-quantized wavefunctions from 
the second-quantized ones, we must "uncompute" the second- 
quantized states. This must be accomplished using only ma- 
nipulations on the register containing the first-quantized wave- 
functions |<^j)i s t: if we can regenerate the input state from 
the wavefunctions, the input register can be reset to |0)2 n d 
as desired. Given the one-to-one correspondence between a 
second quantization state and the corresponding first quan- 
tization wavefunctions, regenerating the input state amounts 
to the problem of identifying the wavefunctions \(j>i)x s t given 
only the information contained in their first-quantized repre- 
sentation. 

For non-degenerate eigenstates, each state \(j>i)i s t can be 
uniquely identified using its energy, which can be obtained 
through the phase estimation procedure ll20l[3ll l32tl . In gen- 
eral, given a unitary U and its eigenstate \ip), the phase es- 
timation algorithm finds the eigenvalue of \ip). Specifically, 
since U\ij}) = e 27 " e |V>),wehave[/ fc |V>) = e 27rife <V). By con- 
trolled applications of the powers of U to controlled on 
the state 5^ EtV |*>. one gets ^ Y,T=o I*) \*) = 
2^/2 T^k=o e 2mke \k) An efficient quantum Fourier 

transform on the control qubits will now yield the first q dig- 
its of the binary expansion of 9. If we choose U such that 
U\4>i)i s t — e 27rlEi \<fii)ist, we can use phase estimation with 
enough control qubits to obtain an approximation of the ener- 
gies Ei. In particular, the natural choice U = e~ lHt , where 
H is the Hartree-Fock Hamiltonian (O, supplies the appropri- 
ate unitary for a suitable choice of the time t. Note that U 
can be simulated efficiently because H is a sum of Fock op- 
erators which are efficiently simulatable by assumption. The 
energy eigenvalues are stored in an additional register con- 
taining enough qubits to provide precision that distinguishes 
between nearby energies. 

In the case that the spectrum of H is degenerate, properties 
other than the energy of the states need to be used to distin- 
guish them. If the degeneracy is caused by a symmetry of the 
Hamiltonian, the elements of the symmetry group can be used 
for this discrimination, as we outline in Sec. Ill El If the de- 
generacies are accidental, other techniques are required, and 
we give some suggestions in Sec. Ill Fl In addition, the tech- 
niques in Sec. Ill Fl can be used for distinguishing eigenvalues 
that are exponentially close together and therefore cannot be 
distinguished efficiently by phase estimation. 

Phase estimation using both U to find energy eigenvalues 
and appropriate symmetry operations to distinguish degener- 
ate states will provide us with a unique combination of eigen- 



values for each state in the superposition. These eigenvalues 
can then be used (for example in conjunction with a look- 
up table) to uniquely identify the wavefunction and subtract 
1 from the corresponding occupation number vector of the 
second-quantization state. Because this is done in superpo- 
sition for every single-particle wavefunction, the input state is 
converted to |0)2nd- 

This accomplishes the total transformation 




which is a separable state. The antisymmetrization step can 
now proceed in superposition as usual, resulting in the final 
state l^ist = J2i a i\^i)ist, as desired. This completes the 
state-preparation algorithm for a given superposition of multi- 
particle states. 

E. Resolving degeneracies caused by symmetry 

The procedure in Sec. IllDl assumes that it is possible to dis- 
tinguish eigenstates based on their energy. If there are degen- 
erate states, additional operations are required to distinguish 
them. Degeneracies in quantum states usually arise as a result 
of symmetry — degeneracies that do not are called "acciden- 
tal" and we treat them separately in Sec. Ill Fl For symmetry- 
caused degeneracy, distinguishing degenerate states requires 
an understanding of how they transform under the symmetry 
operations of the system. All of the wavefunctions \4>i)\ s t 
are eigenstates of each symmetry operation within the point 
group, but degenerate wavefunctions will always have differ- 
ent eigenvalues for at least one of the operations. Phase esti- 
mation can still be used to obtain a unique set of eigenvalues, 
but in addition to finding the energies, we can distinguish the 
wavefunctions by symmetry. By applying phase estimation 
using an appropriate symmetry operation as the unitary opera- 
tor, we obtain additional eigenvalues to distinguish degenerate 
states. 

Because there are only a limited number of symmetries 
that are possible in physical systems, it will rarely be nec- 
essary to use more than a few readout qubits to retrieve all 
the distinguishing eigenvalues. With the exception of sys- 
tems with spherical, cubic, or icosahedral symmetry, which 
we treat below, all systems in three-dimensional space have 
a symmetry point group all of whose irreducible representa- 
tions are one- or two-dimensional 113 311 . Wavefunctions trans- 
forming as the one-dimensional irreducible representations 
are non-degenerate, while the ones transforming as the two- 
dimensional irreducible representations come in degenerate 
pairs. Distinguishing them, therefore, requires the determi- 
nation of only one symmetry eigenvalue which is different for 
the two wavefunctions. 

This is most easily done in the case of point groups C nv , 
Coov, D n , D n h, Dooh, and D n d, all of which contain a Ci 
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axis or a reflection plane that has character zero in all of the 
two-dimensional irreducible representations. In this case, one 
of the two degenerate wavefunctions is invariant under the re- 
flection or C2 rotation, while the other acquires a phase of 
— 1. To distinguish them, one would use the reflection or the 
C2 rotation as the unitary of phase estimation with one read- 
out qubit (note that these operations are easy to implement, 
being simple linear transformations). The readout qubit, ini- 
tialized in the state (|0) + |1)) /v2, would, under the action 
of the symmetry operation, be converted to (|0) ± |1)) /\/2, 
depending on the acquired phase. A Hadamard gate would 
then return |0) or |1), perfectly discriminating between the 
two eigenfunctions. 

Symmetry groups C n , C n h, and S^n have, strictly speak- 
ing, only one-dimensional irreducible representations. How- 
ever, there are pairs of representations that are complex con- 
jugates of each other, meaning that the corresponding energy 
levels are degenerate due to time-reversal symmetry. These 
pairs of conjugate representations are called "separably de- 
generate 13411 ." and the corresponding wavefunctions can be 
distinguished using the principal symmetry axis C n (or Ssn 
in the $2n groups). In eacn case, under the action of C n , one 
of the wavefunctions acquires a phase ui and the other oj*, 
where lj = e 2 ™/" (there are also cases where the pairs ac- 
quire phases such as — u> and — lj*, lj 2 and (w 2 ) , and so on, 
but these do not change the procedure outlined here). Phase 
estimation can, as usual, measure this phase up to a certain 
precision. However, since 1 jn usually does not have a finite 
binary expansion, there will be an associated error in the phase 
estimation. This can be reduced below an arbitrary threshold 
by the addition of more readout qubits, as discussed in Sec. 
IIVI This is especially true since real physical systems almost 
never have C n axes with n > 8, meaning that only several 
qubits will be required for readout. 

The cubic and icosahedral groups, T, Th, Td, O, Oh, I, and 
Ih, all have three-dimensional irreducible representations (and 
I and 1^ also have four- and five-dimensional ones). Fortu- 
nately, there are plenty of reflection planes and C2 axes which 
can be used for discrimination just as was done in the simpler 
groups above. Distinguishing three or four degenerate states 
requires two symmetry eigenvalue comparisons (and three in 
the case of five-fold degeneracy). Consequently, two readout 
qubits are required in these cases, one for each comparison (or 
three qubits in the five-fold degenerate case). 

Degenerate states of spherically symmetric systems, such 
as atoms, can be distinguished by energy and by their angular 
momentum quantum numbers £ and nig. The maximally sym- 
metric case is the 1 jr potential, where the conservation of the 
Laplace-Runge-Lenz vector implies that all states with equal 
principal quantum number n are degenerate. If our basis con- 
tains states with n < n max , we would require 0(log 2 n max ) 
qubits for the discrimination of the angular momentum states 
(that is, 0(log 2 n max ) qubits each for I and mi). While cir- 
cumstances where one encounters states of extremely high an- 
gular momentum are rare, we can see that the discrimination 
can be performed efficiently. The phase estimation in this case 
would use discrete rotations as its unitary operator. A similar 
approach was suggested by Zalka for the related problem of 



implementing unitary representations of SU(2) j35ll . 

F. Resolving accidental degeneracies and exponentially close 
eigenstates 

In Sec. Ill El we outlined a procedure for distinguishing 
states that are degenerate because of symmetry. However, the 
eigenstates might also be accidentally degenerate or exponen- 
tially close in energy so that they cannot be efficiently distin- 
guished by phase estimation. In those cases, it is not possible 
to distinguish between the (near-)degenerate states using the 
symmetry-based procedure. 

One way around these problems is to transform to another 
basis where the (near-)degeneracy does not arise. A way 
of accomplishing this is to use a perturbed Fock operator 
F' = F + V, where V is a small, efficiently simulatable 
perturbation that breaks the (near-)degeneracies. In a finite 
basis, V must also be small to ensure that the new basis can 
adequately describe the target state. The new eigenfunctions 
are obtained from the old using perturbation theory, as are the 
new coefficients of the state that we wish to prepare. This 
change of basis can be done efficiently on a classical com- 
puter, before proceeding as normal with the state preparation 
algorithm. For the purposes of phase estimation, the new Fock 
operator can be efficiently simulated by operator splitting be- 
cause both F and V are efficiently simulatable. 

A drawback of this procedure is that the perturbation may 
destroy certain desirable symmetries of the system. In some 
cases, this can be avoided if we choose V cx \4>i) (4>i\, where 
\4>i) is one of the (near-)degenerate eigenstates. In that case, 
F' and F would have the same eigenstates and no change of 
basis would be needed. Of course, it is possible that V in this 
form is not efficiently simulatable, in which case this scheme 
would not be efficient. 



G. Mixed states 

The previous sections outline the procedure for preparing 
general pure states, which in the chosen basis read 

i«>u t = E* (7) 

i 

From now on, we drop the subscript 1st for clarity. We now 
wish to prepare a mixed state with density operator 

i 

where j^,-) are arbitrary pure states of the form (0 and the 
probabilities pi add up to 1 . This scheme could be used for the 
preparation of thermal states, in which case one would choose 
to be the Hamiltonian eigenstates and = e~" Ei /Z, 
where j3 — l/k^T and Z is the partition function. Our ap- 
proach to the thermalization problem is therefore different 
from that of Terhal and DiVincenzo, who prepare thermal 
states by simulating an external bath 13611 . 
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We assume that each \^>i) can be efficiently specified using 
some specification |&) (for example, |^) is the ^th eigen- 
state of the Hamiltonian). We begin by preparing the state 
Si y/Pi\^i)- This can be done using the procedure in Sec. 
Ill Al if we order the £j's so that they may be thought of as 
a function on a one-dimensional grid. We then run the entire 
state-preparation algorithm in superposition, preparing the ap- 
propriate conditional on the value of the This yields 
the state 

|3> = £^l£i>l*i>, 

i 

the density operator of which is 

Tracing out the specification register, we get the desired den- 
sity operator 

P = Tr 5 p H 

i 

In practical terms, tracing out the specification register 
amounts to doing nothing at all. That is, each is entangled 
to a different meaning that the l^^'s evolve separately 
under time evolution, as they would if they were independent 
members of an ensemble. 

One can notice that density operators diagonal in the 
basis can be prepared more directly. In the previous Sec. 
Ill Dl we had to "disentangle" the first- and second-quantized 
states. If we had instead simply traced out the input register, 
we would have obtained a mixed state diagonal in the 
basis. 



III. MANY TYPES OF PARTICLES 

In Sec. HU we outlined an algorithm for the preparation of 
arbitrary many-particle states (pure or mixed) of a system of 
identical particles. However, one often wants to consider sys- 
tems of more than one type of particle, or treat particles of 
the same kind, but separated in space, as different (the latter 
approach might be useful, for example, in computing electron 
transfer matrix elements for large molecules) l37ll . We con- 
sider the case of two types of particles, with the generalization 
to more types being clear. 

One wants to prepare an arbitrary two-particle state 

\e)=J2<Xi,j\®A,i)\$B,j), 

i,j 

where |3>A,t) is a many -particle eigenstate of particles of type 
A, and \$B,j) is an eigenstate of particles of type B. Each 
element |3>a,i) \®B,j)oi this superposition is easily created 



by preparing the appropriate state in separate registers as was 
done in Sec. Ill CI Creating |0) itself can be done in analogy 
to the preparation of superpositions in Sec 111 Dl We start by 
efficiently specifying |0) using occupation number vectors of 
the |3>A,i) an d the |^s,i)> namely 

^ ai,j\nA,i)2nd\nB,jhnd ® 1st |0s> l st . 

We then complete the state preparation, in superposition, as 
we did in Sec. IUDI treating each register separately. Doing 
so produces |0). 

There are many circumstances in which the ability to pre- 
pare states such as these would be valuable. For instance, in 
chemical dynamics it is necessary to treat the nuclei and the 
electrons separately. If we restricted our state preparation to 
simple product states such as |3>A,i) \®B,j)> we would get a 
state in the Born-Oppenheimer approximation, which is often 
a good approximation to the initial states of reactants partic- 
ipating in chemical reactions. However, as the procedure for 
preparing 1 0) shows, quantum computers could just as easily 
prepare non-Born-Oppenheimer states in which there is cor- 
relation between electronic and nuclear degrees of freedom. 

Many -particle mixed states can likewise be prepared by fol- 
lowing the procedure in Sec. Ill Gl separately for each type of 
particle. 

IV. ERRORS AND THE COMPUTATIONAL COST 

For the state preparation algorithm to be considered effi- 
cient, the time it takes to execute it must scale as a polyno- 
mial in the sizes of the input. More precisely, it should scale 
as a polynomial in I, the number of qubits used to store the 
wavefunction and m, the number of occupied single-particle 
orbitals, which is the best descriptor of the total size of the 
system. 

In this section, we first show that pre-existing errors are am- 
plified at most linearly by subsequent steps of the algorithm. 
We then use this fact to obtain the total computational cost of 
preparing an arbitrary quantum state. 

A. Errors 

Assuming that the quantum gates are executed perfectly — 
or that the gate errors are corrected using efficient error cor- 
rection algorithms — there are five sources of error in the state 
preparation algorithm: 

1. Preparation of single-particle eigenstates. Zalka's 
method that we adopt in Sec. Ill Al requires evaluation of the 
integrals (Q]). We have addressed the computational cost of 
integral evaluation in Sec. Ill 131 where we show that the pro- 
cedure can be accomplished in time polynomial in ej 1 if the 
wavefunction's indefinite integral is known or, more gener- 
ally, if the wavefunction is bounded. The resulting error in the 
prepared single-particle eigenstate is < lei/ 2. 

2. Assembly of many-particle eigenstates. Many -particle 
eigenstates (01 inherit the errors present in the single-particle 
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eigenstates \fa) that are used to assemble them. Supposing 
that the prepared states approximate the true states \<fii) 

with error = 1 — (fa fa} , then the prepared Hartree 

product 4>i, . . . fa_ > suffers an error 



e$ = 1 - 

= l-JJ(l- £ ^,) 

i=l 

m 

< ^ £j>,ji < me^ < mlei/2, 

i=l 

where = maxe^ ;i . Since the total error grows as a poly- 
nomial in both the single-state error and the number of oc- 
cupied states, the assembly of Hartree products amplifies the 
pre-existing errors only linearly in m. The remaining step, the 
antisymmetrization of the Hartree product, does not introduce 
additional errors. 

3. Preparation of superpositions. The parallel state- 
preparation that is used to perform the transformation ((6) does 
not introduce any additional errors with the exception of the 
possible failures of phase estimation, discussed below. Nev- 
ertheless, we should see how pre-existing errors propagate 



through this step. If the prepared state is 



4\ 



we see that it suffers an error with respect to the target state 



4\ 



where e$ 



4>, 



i 

= l-^loil^l-e^O 

i 

< £$ < mlei/2, 
- maxe^ i and where we have assumed that 



j / — for i ^ j. In other words, the error in 

is limited by the error of its components. 

4. Discrimination of states in a superposition. The prepa- 
ration of superpositions described in Sees. Ill DHII Fl andUIIlre- 
lies on phase estimation as a means of distinguishing states. 
Since the eigenenergies will rarely have finite binary expan- 
sions, there will be errors introduced at this step. If two phases 
differ at the nth bit and we perform phase estimation with 
q = n + p qubits, the probability of an incorrect identifi- 
cation is 1/2(2 P — 2), meaning that the success probability 
will be 1 — £pe provided we implement phase estimation with 
p = [log (2 + l/2e PE )l additional qubits HI. The addi- 
tional overhead, logarithmic in £pg, does not compromise the 
efficiency. The same arguments apply to the phase estimation 
of eigenvalues of C n belonging to states in separably degen- 
erate irreducible representations of groups C n , C n h, and S^n 
(see Sec. Ill Eb . The symmetry eigenvalues that are useful for 
states in the other point groups are always ±1, and can be per- 
fectly resolved using phase estimation with a single readout 
qubit. 

In addition, failures of state discrimination based on phase 
estimation can be detected after the state preparation is com- 



plete. The second-quantized register, which should be uncom- 
puted during the procedure, should be measured at the end. 
If |0 . . . ) is observed, phase estimation will have succeeded. 
Otherwise, a misidentification will have occurred, and the pro- 
cedure ought to be repeated. This simple, classical error cor- 
rection introduces only a constant overhead and ensures that 
phase estimation does not contribute to the error in the final 
prepared state. 

5. Assembly of mixed states. In the notation of Sec. Ill Gl 



if the prepared states 



with an error e^j = 1 — 



^i) approximate the true states 



and assuming per- 



fect preparation of the state J2i \/Pi l&}> tne final prepared 
mixed state will be p = J^i Pi 



*i > < * 



. If we assume that 



=0 for i 7^ j, then p suffers an error Hill 



£ P = 1 - Tr V PPV P 



1/2 



1/2 



i 

< < mlei/2, 

where — max£$.,. That is, the assembly of mixed states 
does not magnify the pre-existing errors. 

Overall, we see that errors introduced in any stage of the 
state preparation algorithm are not amplified more than poly- 
nomially by subsequent stages. The final error in the prepared 
state is e = e p < mlei/2, meaning that the error scales lin- 
early with the size of the system m and the error of the in- 
tegration procedure, as well as logarithmically with the grid 
size 2 l . 



B. Computational cost 

There are three time-consuming steps in the state prepa- 
ration algorithm. The first is the evaluation of the integrals 
(Q3 and the resulting single-qubit rotations, the second is the 
phase-estimation that is used to distinguish states in the super- 
position (see Sec. Ill Dl ). and the final is the antisymmetrization 
procedure described in Sec. Ill CI We characterize the cost of 
each step in turn. 

In the previous section, we have seen that the total error of 
the prepared state will be e < mlej/2. Therefore, if we want 
to ensure a maximum error e, we must choose ei = 2s /ml, 
implying that 0(mZe _1 ) time is required for each integration 
(see Sec. HI BP . The integration procedure itself is called ml 
times: for each of the m occupied orbitals, / qubits have to be 
rotated correctly. Therefore, the total time required for all the 
qubit rotations is 0(m 2 l 2 e~ 1 ). 

The cost of the phase-estimation procedure that is used to 
distinguish the eigenstates cannot be given precisely because 
we have not made any assumptions about the nature of the 
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Fock operator F other than that it is efficiently simulatable, 
that is, running in time poly(m, M , I, A -1 ) (here, A is the 
precision at which the simulation needs to be run, i.e., it is 
half the gap between the closest two eigenstates, which we 
assumed is not exponentially small). Simulating the entire 
Hartree-Fock Hamiltonian requires the simulation of the Fock 
operator acting separately on each particle, meaning that the 
total simulation requires mpoly(m, M, I, A -1 ) time. In ad- 
dition to this, two quantum Fourier transforms (QFTs) are re- 
quired on the readout register of the phase estimation. If q 
qubits are used for the readout (see Sec. |IV A| 4), the QFTs 
require 0(q 2 ) time. It should be noted that the required q 
is determined only by needed precision in the phase estima- 
tion, and that it does not depend strongly on m, I, or M. 
Therefore, the cost of the QFTs can be treated as essentially 
a constant overhead. Furthermore, there is the cost of look- 
ing up the state's energy in the look-up table; a simple binary 
search requires 0(log 2 M) time per register, for a total cost 
of 0(mlog 2 M). But this, too, is a negligible cost in com- 
parison to mpoly(m, M, I, A -1 ), which we conclude is the 
asymptotic cost of the eigenstate discrimination portion of the 
state preparation algorithm. 

The bottleneck of the antisymmetrization procedure used to 
produce fermionic states (or the symmetrization for bosonic 
ones) is the sort that takes state (0]i to (0. Sorting register 
B by a comparison sort requires f2(mlogm) swaps. These 
swaps must also be performed on each of the corresponding I 
qubits of register A, for a total cost of ft(lm log m). For large 
systems, this expression will be dominated by the scalings of 
the integral evaluation and the phase estimation. 

Based on the foregoing, the total computational cost 
of the state preparation algorithm is 0(m 2 l 2 e^ 1 + 



mpoly(m, M, I, A -1 )) = poly(m, M, I, A -1 ), an ex- 
pression polynomial in all the basic descriptors of the system. 
This allows us to conclude that the algorithm, as described 
above, is efficient. 



V. CONCLUSION 

We have outlined a quantum algorithm for the preparation 
of physically realistic quantum states on a lattice. In particu- 
lar, we have gone beyond previous proposals by describing a 
method for preparing any pure or mixed state of any number of 
particles. This is achieved by using Zalka's method for prepar- 
ing single-particle states and then combining those into many- 
particle states. The assembly of many-particle states requires 
that we be able to distinguish them on a quantum computer, a 
task that we address using phase estimation. We also provided 
symmetry-based solutions for degenerate cases, where phase 
estimation using a single operator is insufficient to distinguish 
the states. Accidentally degenerate states can be distinguished 
by adding a perturbation to the system Hamiltonian. Our algo- 
rithm is efficient, with a run-time of poly(m, M, I, e -1 , A -1 ), 
subject only to the requirements that the wavefunction be 
bounded or that its indefinite integral be known and that the 
Fock operator be efficiently simulatable. 
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